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1.  Introduction 


This  report  presents  a  set  of  functions,  written  in  C++,  that  is  designed  to  work  with  evenly 
spaced  rectangular  surface  grids.  Grids  of  this  type  have  a  variety  of  applications,  including 
representing  terrain  features  and  storing  spatial  probability  information.  Relative  to  other 
methods  of  storing  surface  information,  they  can  be  advantageous  when  calculation  time  is 
critical. 

The  functions  that  are  described  in  this  report  have  been  grouped  into  the  yBilinear  namespace, 
which  is  summarized  at  the  end  of  this  report.  Functions  for  calculating  surface  interpolations, 
line-surface  intersections,  and  surface  gradients  are  included.  The  functions  are  based  on  bilinear 
interpolations  and  have  been  templated  to  allow  for  a  variety  of  pointers  and  containers. 

The  yBilinear  namespace  relies  exclusively  on  standard  C++  operations  and  functions.  However, 
example  code  that  is  included  in  this  report  makes  use  of  the  yRandom  namespace-  for 
generating  pseudorandom  numbers  and  the  yBmp  namespace-  for  creating  images. 


2.  Derivations 


2.1  Evenly  Spaced  Rectangular  Surface  Grids  Defined 

Surface  information  can  be  stored  in  evenly  spaced  rectangular  grids,  as  shown  in  Fig.  1,  where 
each  dot  represents  a  grid  point  with  associated  x,  y,  and  z  coordinates. 

Since  surface  grids  are  assumed  to  be  evenly  spaced,  x,  and  y ,  values  can  be  stored  implicitly: 

x,  =  x0  +  /Ax  for  0  <  /  <  m  ( 1) 


and 

y  j  =  To  +  J'Ay  for  0  ^  j  <  n ,  (2) 

where  x0  and  y0  are  used  to  locate  the  lower-left  comer  of  the  grid,  and  Ax  and  Ay  are  defined 
to  be  the  difference  between  consecutive  x  and  y  values,  respectively. 
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Fig.  1  An  evenly  spaced  rectangular  surface  grid 


Eqs.  5  and  6  can  be  used  to  find  grid  indices  given  spatial  coordinates. 

/ 

i  =  trunc 


^  x  —  x0  ' 


and 


j  =  trunc 


Ax 


f  \ 

y-y  o 


v  Ay  j 

where  the  trunc()  function  rounds  toward  zero  to  the  nearest  integer. 

Eqs.  7  and  8  can  be  used  to  find  maximum  values  for  the  x  and  y  coordinates. 

*max  =*„-!  =X0+(m-l)Ax. 

fmax  =y„-i  =y0+(n-i)Ay. 


(5) 

(6) 


(7) 

(8) 


2 


2.2  Bilinear  Interpolations 

Bilinear  interpolations  can  be  used  to  estimate  surface  values  for  locations  that  are  between  grid 
points. 

To  calculate  a  bilinear  interpolation,  begin  with  the  equation  for  calculating  a  linear 
interpolation,  which  is  easily  derived  from  the  point-slope  equation  of  a  line  (Fig.  2): 


(9) 


b  =  ya-axa 


a 


Fig.  2  A  linear  interpolation 

A  bilinear  interpolation  can  be  calculated  by  performing  3  linear  interpolations.  First,  use  Eq.  9 
to  interpolate  in  the  y  direction  (shown  in  green  in  Fig.  3): 


(10) 


(11) 


Next,  interpolate  in  the  x  direction  (shown  in  red  in  Fig.  3): 


(12) 
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Fig.  3  A  bilinear  interpolation 

Eqs.  2  and  4  can  be  used  to  rewrite  Eqs.  10  and  1 1  in  terms  of  y0  and  Ay  : 


z„  = 


'y-y0 

v  Ay 


•/ W/-I- 


Z  n  = 


y-y° 


Ay 


j  !  i+Uj  ‘ 


Similarly,  Eqs.  1  and  3  can  be  used  to  rewrite  Eq.  12  in  terms  of  x0  and  Ax 


z  = 


y  x  —  X, 
v  Ax 


--^(Zfi~Za)+Za- 


2.3  Line-Surface  Intersections  (Cells) 

Suppose  that  a  line,  L  ,  passes  through  the  points  Z0  and  Lx ,  where 

h  =  LoJ  +  LoJ  +  LoJ  and  A  =  LiJ  +  A Jy  +  LiJ  ■ 

Let  p  represent  a  point  that  lies  on  L  ,  where 
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p  -xx  +  yy  +  zz . 


(17) 


L0  and  Lx  can  be  used  to  construct  a  parametric  equation  for  p  as  a  function  of  t : 

p  =  (Ll-L0)t  +  L0.  (18) 

The  parameter  t  represents  the  scaled  distance  from  Z0  to  Lx  along  L  .  Thus,  if  t  =  0 ,  p  is 
located  at  L0 .  If  t  =  1 ,  p  is  located  at  Lx . 


From  Eqs.  17  and  18, 

x  =  C0t  +  L0x  ,  y  =  Cxt  +  L0y  ,  and  z  =  C2t  +  L0z,  (19) 

where 

C0  =  Ll  x  —  L0  x  ,  Cj  =  Lly  —  L0  y  ,  and  C2  =  LX  z  —  L0  z .  (20) 

Substituting  Eqs.  10  and  1 1  into  Eq.  12  and  making  use  of  Eqs.  3  and  4, 

z(x,  y)  =  [(x-  xM  )( y  -  yj+l  )ztj -(x- xt )( y  -  yj+1  )zi+lj  ^ 

-(x- XM X v - yj)zij+l  +  (x- Xj X v -  Vj )zMJ+l ] /  AxAv . 

Substituting  Eq.  19  into  Eq.  21  results  in  an  equation  that  can  be  used  to  find  t  at  the  point  where 
L  intersects  the  surface: 

C2t  +  L0z  =  [(C0t  +  L0x  —  xM  ){Cxt  +  L0y  —  y  j+l  )ziy 
~ (Q  +  L0x  - x,)(Q  +  L0y  - yj+1)zMJ 
~ (Q  +  L0x  - xM )(Cxt  +  L0y  - yj)zi  j+l 
+  (C0t  +  L0x  -  x,  )(C\l  +  L0y  -  yj  )zMJ+l  ]  /  Ax  Ay . 

Rearranging  terms, 

(C2t  +  L0  z ) AxAy  =  (C0t-  (xi+l  -  L0x )){Cxt  -  (yj+x  -  L0  y ))zt  j 
- (C0/  - (x,  - L0x))(Cxt  - (y j+l  -  L0y))zj+l  j 
-  (C0f  -  (x,.+1  -  L0  x))(Cxt  -  (yj  - L0y))zi  j+X 
+  ( CQt  -  (x;.  -  L0  x))(Cxt  -  (yj  -  L0y))zj+X  j+1 . 

Another  set  of  constants  can  be  used  to  make  Eq.  23  a  bit  more  manageable: 

C3  =  A+ 1  —  ^o,x  5  C4  =  y j+x  —  L0  y  ,  C5  =  x,  —  L0  x  ,C6  =  y y-  —  T0  ^ . 


(22) 


(23) 


(24) 


Thus, 
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(C2t  +  L0z)AxAy  =  ( C0t  -  C3)(Q  -  C4)ztJ 

~  (C0t  -  C5  ){Cxt  -  C4  )zMj 
-(Q-C3)(Q-C6)z,.+1 
+  (C0t-C5)(C1t-C6)zMJ+l. 

Expanding  quadratics, 

(C2t  +  L0z)AxAy  =  (C0Q2 ~(ClC3+C0C4)t  +  C3C4)zlJ 

-(QQ2  - (ClC5  +  C0C4)t  +  C5C4)zi+lj 
-(QQ2  -  (ClC3  +  C0C6)t  +  C2C6)zij+l 
+  (C0Cjt2  -  (C,C5  +  C0C6)t  +  C5C6)zi+l  j+1 . 

Regrouping  terms, 

CoQCz,.- 

Zi+\,j  Zi,j+\  +Z;+lJ+l)^ 

+[-(C1C3  +  C0C4)z, .  +  (CjC5  +  C0C4)z,+u 
+  (CjC3  +  C0C6)z,  ,+1  -  (C3C5  +  C0C6)zi+u+1  -  C2AxAv]t 
_i_^,3^'4'z/,7  C5C4zi+iJ  —  C3C6z(.^.+1  +  C5C6zM  J+l  —  L0zAxAy  =  0 . 


Note  that  Eq.  27  is  a  quadratic  in  standard  form: 

at2  +  bt  +  c  =  0 , 

where 
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a  =  C0C,  (z,.j  -  zMJ  -  z.  y+1  +  zi+u+1 ), 

-(CtC3  +  C0C4)z, .  +  (CjC5  +  C0C4)z,.+u 
+  (C,C3  +  C0C6)z,,+1  -  (C,C5  +  C0C6)z,+u+1  -  CAvAv  , 


and 


c  —  C3CAzij  C5CAzMj  C3C6zij+l  +  C5C6zMj+1  L0 


Ax  Ay . 


Finally,  solving  for  t , 


5  |  5  e  9? 
undefined 
c 
b 

-  b±  ^Jb2  -  4oc 
2  a 


for  a  =  0,  b  =  0,  and  c  =  0 
for  a  =  0,  b  =  0,  and  c  ^  0 

for  a  =  0  and  b  ^  0 

for  o  0 


(25) 


(26) 


(27) 


(28) 

(29) 

(30) 


(31) 


(32) 
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For  the  a  ^  0  case,  if  the  discriminant  is  less  than  zero,  then  no  real  solutions  exist,  implying  that 
L  does  not  intersect  the  surface. 

2.4  Line-Surface  Intersections  (Grids) 

Testing  for  intersection  between  a  line  segment,  L  ,  and  an  evenly  spaced  rectangular  surface 
grid  involves  checking  for  intersection  between  one  or  more  simple  cell  surfaces,  each  defined 
by  Eq.  21.  Figure  4  shows  the  projection  of  L  onto  the  x  —  y  plane,  with  the  cells  that  L  might 

intersect  highlighted  in  yellow. 


Fig.  4  Possible  cell  intersections  for  line  segment  L  intersecting  a  surface 


The  following  algorithm  can  be  used  to  find  indices  that  are  associated  with  cells  that  line 
segment  L  might  intersect. 

1.  Initial  Calculations: 

a.  Define  i  and  j  to  be  the  x  and  y  indices,  respectively,  of  the  current  cell  that  is 
being  checked.  Use  Eqs.  5  and  6,  along  with  values  for  L0  x  and  L0  v ,  to  set  initial 
values  for  i  and  j  .  Restrict  initial  values  of  i  |  0  <  i  <  m  - 1  and  j  |  0  <  j  <  n  - 1 . 


7 


b.  Define  /  and  J  to  be  the  x  and  y  indices,  respectively,  that  are  associated  with  the 
end  of  the  line  segment.  Use  Eqs.  5  and  6,  along  with  values  for  LX  x  and  L  ,  to  set  / 
and  J  .  Restrict  values  of  /  I  0  <  /  <  m  - 1  and  J\0<J<n-\. 


c.  Define  A i  and  A /  to  be  directional  quantities  that  can  be  used  to  increment  i  and  j  . 
Use  Eqs.  33  and  34  to  set  A i  and  A j  : 


A  i 


I 

-1 


4/  = 


l 

-l 


for  i  <  I 

otherwise 

(33) 

for  j  <  J 

otherwise 

(34) 

2.  Outer  Loop: 

a.  Define  i,  to  be  the  index  associated  with  the  location  where  the  projection  of  L  onto 
the  x—y  plane  crosses  a  constant-v  reference  line  (e.g.,  the  location  of  the  green  dot  in 
Fig.  4). 

b.  If  j  =  J  ,  set  it  =  I . 

Otherwise,  use  Eqs.  35  and  36  to  calculate  yt  then  xt .  Use  xt ,  along  with  Eq.  5,  to 
calculate  i, .  Restrict  the  value  of  it  \  0  <  it  <  m  - 1 . 


j  (j  + 1)4 y  +  y0  for  A j  =  1 
[  jAy  +  y0  otherwise 


xt  = 


y,  -  A 


O.v 


Ll,y  L0  ,y 


(A,x  Lq  x)+  L0  x 


(35) 

(36) 


3.  Inner  Loop: 

a.  Store  i  and  j  . 

b.  If  Ai(i  -  it )  <  0  ,  increment  i  by  A i  then  repeat  step  3a.  Otherwise,  continue  to  step  4. 

4.  End  Condition: 

a.  Increment  j  by  A j  and  i  by  -  A i . 

b.  If  A j(J  -  /)  >  0  ,  then  jump  back  to  step  2b.  Otherwise,  the  algorithm  is  complete.  The 
stored  values  for  i  and  j  represent  the  complete  set  of  cells  that  line  segment  L  may 
intersect. 
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2.5  Surface  Gradients  and  Normal  Vectors 

To  avoid  confusion  between  the  dependent  variable  z(x,y)  and  the  independent  variable  z  , 
define  </>(x,y)  =  z(x,y) .  Then,  from  Eq.  21, 


<p(*>  y)  =  [(*  -  *,-+i  )(y  -  yj+i  Kj  -  O  -  a  )(y  -  yJ+i 


M  ,j 


-(x-  xM  )(y  ~  Vj  )ziJ+1  +  (x-xi)(y-  Vj  )zMJ+1  ]  /  AxAv . 
From  the  definition  of  the  del  operator, 

,  df  „  df  ,  df  , 

Vf(x,y,z)=—x  +—y  +  —z 
ox  oy  oz 

=>V<p  =  [(y-  yj+l  )(zij  -  zMj  )  /  Ax  Ay  -  ( y  -  v;.  )(zij+1  -  zi+1J+l )  /  Ax  Ay]  x 
+  [(x~  xM  )(zu  -  ziJ+1 )  /  Ax  Ay -(x- x,  )(zi+lJ  -  zMJ+1 )  /  Ax  Ay]  y . 


(37) 


(38) 

(39) 


Thus,  Eq.  39  is  the  gradient  of  the  surface  defined  by  Eq.  21.  Note  that  is  not  of  unit  length. 
A  number  of  useful  surface  properties  can  easily  be  calculated  based  on  the  gradient;- 

The  x  —  y  components  of  the  direction  of  steepest  ascent  are  given  by  (V  (j))x  and  (V  (j))  v  from 
Eq.  39. 

The  slope  in  the  direction  of  steepest  ascent  is  given  by 


V^=J(V^)x2+(V^)v2 


and  the  unit-nonnal  vector  n  is  given  by 

-WU-(V0,j>  +  £ 


n  = 


(v^)x2+(m,2+ 1 


(40) 


(41) 


Recall  that  unit-normal  surface  vectors  are  not  unique.  Specifically,  if  h  is  a  unit-normal  surface 
vector,  then  -  n  is  also  a  unit-normal  surface  vector,  n  has  been  chosen  to  have  a  positive  z 
component. 

2.6  Scaled  Coordinates 

Eqs.  42  and  43  provide  a  means  to  convert  to  and  from  scaled  coordinates. 


sr  = 


^and 


\v  7  Ay 

x  =  x0  +  sxAa  and  y  =  y0  +  s  Ay . 


(42) 

(43) 
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The  benefit  of  working  with  scaled  coordinates  is  that  they  simplify  calculations  by  introducing 
unit  intervals  and  zero  offsets.  Thus,  they  allow  for  making  the  following  substitutions  when 
working  with  equations  prior  to  Eq.  42: 

x  -»  ,  y  -»  sy ,  Ar  ->  1 ,  by  -»  1 ,  x0  -»  0 ,  y0  ->  0 ,  xt  ->  i ,  and  y}  ->  j .  (44) 

The  use  of  scaled  coordinates  in  the  code  that  follows  has  2  advantages:  they  make  the  code  run 
slightly  faster,  and  they  reduce  the  required  number  of  user-supplied  parameters. 


3.  Calculating  Indices:  The  SafelndexQ  Function 


The  Safelndex()  function  uses  Eq.  5  (or  6),  as  well  as  the  substitutions  from  Eq.  44,  to  calculate 
the  grid  index  that  is  associated  with  a  user-supplied  scaled  coordinate.  The  grid  index  is  then 
modified  to  be  compatible  with  other  functions  described  in  this  report: 


0 


k  =  < 


m  —  2 
trunc(.s) 


for  s  <  0 
for  s  >  m  -  2 . 
otherwise 


(45) 


3.1  SafelndexQ  Code 


inline  int  Safelndex(//<=====CALCULATE  A  SAFE  INDEX  AT  A  GIVEN  SCALED  LOCATION 

int  m, //< . NUMBER  OF  GRID  INDICES  (m  SHOULD  BE  AT  LEAST  2) 

double  s ){//< . A  GRID  LOCATION  IN  SCALED  COORDINATES 


return  s<0?0: s>m-2?m-2 : int(s) ; 

}//  rj  rsj  rsjrsj  YAGENAUT@GMAIL.COM - LAST~UPDATED~16AUG2014 - 

3.2  Safelndex()  Parameters 

m  m  specifies  the  array  size  that  is  associated  with  the  scaled  distance  s.  Typically, 

m  refers  to  m  when  s  =sx  and  n  when  s  =s  .  m  should  be  greater  than  1. 

s  s  specifies  a  scaled  distance,  typically  either  sx  or  sy  from  Eq.  42. 

3.3  Safelndex()  Return  Value 

The  SafelndexQ  function  returns  the  index  k  from  Eq.  45. 
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4.  Performing  Interpolations:  The  Interpolate!)  Function 


The  Interpolate()  function  uses  Eqs.  13,  14,  and  15,  as  well  as  the  substitutions  from  Eq.  44,  to 
perform  bilinear  interpolations.  The  equations  have  been  coded  in  a  slightly  more  compact  fonn: 

Z«  =  (Sy  -  j\zij+ 1  -  ZiJ  )+  Zi,j  '  (46) 

Z  =  -  0{k  -  j\Zi+l,j+\  -  ZM ,j)+  ZMJ  -  Za  }+  Z«  '  (47) 

4.1  Interpolate()  Code 

template<class  T>double  Interpolate(//<=================BILINEAR  INTERPOLATION 

const  T&Z,  //< . POINTER  TO  Z  WHERE  Z [ i]  [  j  ]=Z(X[i]  ,Y[  j  ]  ) 

int  ij  int  j,//< . GRID  INDICES  (0<=i<m-l  &  0<=j<n-l) 

double  sx, double  sy){//< . A  GRID  LOCATION  IN  SCALED  COORDINATES 

double  za=(sy- j)*(Z[i] [ j+1] -Z[i] [ j ] )+Z[i] [j ] ; 

return  (sx-i)*((sy-j)*(Z[i+l] [ j+l]-Z[i+l] [ j])+Z[i+l] [ j]-za)+za; 

}//  YAGENAUT@GMAIL.COM - LAST-UPDATED-16AUG2014 - 


4.2  Interpolate()  Parameters 

Z  Z  points  to  a  2-index  array  of  values  that  represents  z(x,y) ,  where 

Z[i][j]  =  z(Wy)- 

i  i  specifies  i ,  an  index  that  is  associated  with  the  x  axis.  Values  for  i  cannot  be 

less  than  zero  or  greater  than  m-2,  where  in  is  the  first-index  size  of  Z. 

j  j  specifies  j  ,  an  index  that  is  associated  with  the  y  axis.  Values  for  j  cannot  be 

less  than  zero  or  greater  than  n  -  2  ,  where  n  is  the  second-index  size  of  Z. 

sx  sx  specifies  sx  ,  a  scaled  distance  from  Eq.  42. 

sy  sy  specifies  s  ,  a  scaled  distance  from  Eq.  42. 

4.3  Interpolate!)  Return  Value 

The  Interpolate!)  function  returns  z  from  Eq.  47. 

4.4  Interpolate!)  Simple  Example 

The  following  example  begins  by  defining  the  surface  that  is  shown  in  Fig.  5.  The  Interpolate!) 
function  is  then  used  to  estimate  z  at  x  =  1 .5  and  y  =  0.5 . 
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#include  <cstdio>// . printf() 

#include  "y  bilinear .  h"// . yBilinear 

int  main( ){//<=================A  SIMPLE  EXAMPLE  USING  THE  Interpolate( )  FUNCTION 

const  int  m=5,n=3; 

double  Z[m][n]={l,l,0  ,  1,0,0  ,  0,1,0  ,  0,0,0  ,  1,0,1}; 
double  x=1.5,y=.5; 

int  i=yBilinear: :SafeIndex(m,x), j=yBilinear: : Safelndex(n,y) ; 
double  z=yBilinear: :Interpolate(Z,i, j,x,y); 
printf("At  x=%.lf  and  y=%.lf,  z=%.lf\n",x,y,z); 

^  j  j  <v»v/vivfvivY  AG  E  N  AUT  ^)GF1A  I  LAST  r^j  LJ  P  DAT  E  D  ^  1 6  AUG  2  0 1. 4^^° 


OUTPUT: 

At  x=1.5  and  y=0.5,  z=0.5 


4.5  Interpolate()  Image  Example 

The  following  example  begins  by  defining  the  Rainbow()  function,  which  can  be  used  to  map 
numbers  to  colors.  Next,  the  example  defines  the  surface  that  is  shown  in  Fig.  5.  Functions  from 
the  yBmp  namespace,  along  with  the  Interpolate()  function,  are  used  to  create  a  pseudo-color 
image  of  the  surface.  Finally,  the  image  is  written  to  a  BMP  file,  which  is  displayed  in  Fig.  6. 
The  color-to-value  key  is  shown  in  Fig.  7. 
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inline  void  Rainbow(//<========================================RAINBOW  COLOR  MAP 

unsigned  char  C [3] , //< . OUTPUT  COLOR  (CALCULATED) 

double  x,//< . VALUE  FOR  WHICH  A  COLOR  WILL  BE  CALCULATED 

double  min, double  max){//< . MINIMUM  AND  MAXIMUM  SCALED  VALUES 

if(x<min){C[0]=C[l]=C[2]=0;/*&*/return;}// . set  too  small  values  to  black 

if (x>max){C[0]=C[l]=C[2]=255;/*&*/return; }// . set  too  large  values  to  white 

x=(l-(x-min)/(max-min))*8;// . remap  x  to  a  range  of  8  to  0 

C[0]=int( (3<x&&x<5 | |x>7  ?-fabs(x/2-3)+1.5:5<=x&&x<=7?l:0)*255);// . blue 

C[l]=int( (l<x&&x<3 1  j  5<x&&x<7?-fabs(x/2-2)+1.5:3<=x&&x<  =  5?l:0)*255);// _ green 

C[2]=int( (  x<l j  | 3<x&&x<5?-fabs(x/2-l)+l . 5 : 1<=x&&x<=3?1 :0)*255) ; // . red 

}//—  ~~YAG  E  N AUT@GMAI L . COM~~~~~~~~~~~~~~~~~~~~~~~~~LAST~UPDATED~04IUN2014-~~~~~~ 


#include  "y_bmp.h"// . yBmp 

#include  "y  bilinear .  h"// . yBilinear, <cmath>{f abs( ) } 

int  main( ){//<================CREATE  AN  IMAGE  FROM  A  SURFACE  USING  Interpolate^) 

const  int  m=5, n=3,M=1000,N=500; 

double  Z[m] [n]={l,l,0  ,  1,0,0  ,  0,1,0  ,  0,0,0  ,  1,0,1}; 
unsigned  char*I=yBmp: :NewImage(M,N,255); 
for(int  q=0;q<M;++q)for(int  p=0;p<N;++p){ 

double  x=q*(m-l)*l./(M-l) ,y=p*(n-l)*l./(N-l); 

int  i=yBilinear: :SafeIndex(m,x), j=yBilinear: :SafeIndex(n,y); 

double  z=yBilinear: : Interpolate (Z, i, j,x,y); 

Ra inbow (yBmp : : Get Pixel ( I, q, p) , z,0, 1) ; } 
yBmp : :WriteBmpFile( "interpolate . bmp", I); 

^  j  j  /v»va/ivmmY AG  E  N AUT ^)GF1A XL*  LAS  T^U P  DAT  E  D 1 6 AUG  2  0  X 4^^ °° 


Fig.  6  Image  generated  by  example  code  from  Section  4.5. 


Fig.  7  Color-to- value  key 
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5.  Calculating  Line-Cell  Intersections:  The  Celllntersect()  Function 


The  Celllntersect()  function  uses  equations  from  Section  2.3  to  calculate  the  point  where  a  line 
intersects  a  simple  surface  that  is  defined  by  Eq.  21.  To  simplify  the  Celllntersect()  function,  the 
equations  from  Section  2.3  have  to  been  rewritten  in  terms  of  scaled  coordinates: 

Eqs.  20  and  24  can  be  converted  to  scaled  coordinates  using  the  substitutions  given  in  Eq.  44: 

C0  =  L\s*  —  L0s^  ,  C,  =  Lx  ^  —  L0^  ,  C2  =  Llz  —  Lqz  .  (48. 

C2=i  +  l-L0  ,  C4=j  +  1-L0  ,  C5=/-E0  ,  C6=j-LQ  (49; 


The  same  is  true  for  Eqs.  29,  30,  and  3 1 : 

a  =  C0Cl(zij-zMj-zu 


j'+l  , 


b  =  -(CjC3  +  C0C4)Z, .  +  (CtC5  +  C0C4)zMj 

+  (C>C3  +C0C6)Z.  ,+1  - (CjC5  +  C0C6)zj+1J+1  - C2 , 


c  —  C2C4ztj  C5C4zj+l  j.  C3C6zi  j+x  +  C5C6zi+l  j+1  L0  z .  (. 

Scaled  coordinates  at  the  point  of  intersection  can  be  found  by  substituting  t  from  Eq.  32  into 
Eq.  19  (after  converting  to  scaled  coordinates): 

sx  —  C0t  +  L0  ^  ,  sy  =  Cxt  +  E0  ,  z  —  C2t  +  L0z .  (. 


Note  that  line  L  will  intersect  the  surface  at  0,  1,  2,  or  infinitely  many  locations. 


If  L  intersects  the  surface  at  2  locations,  then  the  rules  for  selecting  which  intersection  to  choose 
are  somewhat  complicated:  The  Celllntersect()  function  attempts  to  find  the  first  intersection, 
when  traveling  from  L0  toward  Z, ,  that  is  within  the  bounds  given  by  xt  <  sx  <  xM  and 
y :  <  s  <  yj+l .  If  no  intersection  is  found,  then  the  function  attempts  to  find  the  first  intersection 

when  traveling  from  L0  toward  Lx  (regardless  of  whether  or  not  it  is  within  the  cell’s 
boundaries).  If  an  intersection  still  isn’t  found,  then  the  function  uses  the  first  intersection  when 
traveling  from  L0  away  from  Lx . 
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5.1  Celllntersect()  Code 


template<class  T>int  Celllntersect(//<==================LINE-CELL  INTERSECTION 

const  T&Z,//< . POINTER  TO  Z  WHERE  Z[ i]  [  j  ]  =Z(X[ i] ,  Y[  j  ] ) 

int  i,int  j,//< . GRID  INDICES  (0<=i<m-l  &  0<=j<n-l) 

const  double  L[6],//< . A  SCALED  LINE  {L0SX, L0SY, L0Z, L1SX, L1SY, LIZ} 

double&sx,double&sy,double&z,//< - SCALED  INTERSECTION  POINT  (CALCULATED) 

double  epsilon=lE -9 ){//<- . -VALUE  FOR  DIVISION-BY-ZERO  CHECK 

double  C0=L[3]-L[0],C1=L[4]-L[1],C2=L[5]-L[2]; 
double  C5=i-L[0],C6=j-L[l],C3=C5+l,C4=C6+l; 
double  a=C0*Cl*(Z[i] [ j ] -Z[i+1] [ j ] -Z[i] [ j+l]+Z[i+l] [ j+1] ) , 
b=-(Cl*C3+C0*C4)*Z[i] [ j ]+(Cl*C5+C0*C4)*Z[i+l] [j] 

+(Cl*C3+C0*C6)*Z[i] [ j+l]-(Cl*C5+C0*C6)*Z[i+l] [ j+l]-C2, 
c=C3*C4*Z[i] [ j ] -C5*C4*Z[i+l] [ j ] -C3*C6*Z[i] [ j+l]+C5*C6*Z[i+l] [ j+l]-L[2]; 
double  tjD; 

if (tabs (a )<epsilon&&fabs(b)<epsilon) return  tabs (c)<epsilon?- 1:0; 
if (fabs(a)<epsilon)t=-c/b, sx=C0*t+L[0] , sy=Cl*t+L[l] , z=C2*t+L[2] ; 
else{ 

if ( (D=b*b-4*a*c)<0)return  0;//..no  intersection,  <sxJsy,sz>  not  calculated 
t=( -b-sqrt (D) )/2/a , sx=C0*t+L[0] , sy=Cl*t+L[l] , z=C2*t+L[2] ; 
double  tb=( -b+sqrt(D) )/2/a, sxb=C0*tb+L[0] , syb=Cl*tb+L[l] ,zb=C2*tb+L[2] ; 
if (t>=0&&t<=l&&tb>=0&&tb<=l){ 

bool  ra=sx<i | | sx>i+l | | sy< j | | sy> j+l?l : 0; 
bool  rb=sxb<i j  | sxb>i+l | | syb< j | | syb> j+l?l : 0; 
if (ra&&rb&&t>tb){sx=sxb, sy=syb,z=zb; return  1; } 
else  if ( ra | | ! rb&&t>tb)t=tb, sx=sxb, sy=syb, z=zb; } 
else  if (t>tb&&tb>0 1 | t<=0&&t<tb)t=tb, sx=sxb, sy=syb, z=zb; } 

if  (sx<i  |  |  sx>i+l  |  |  sy<  j  |  |  sy>j+l)return  1;// . out  of  bounds 

return  t<0?2 : (t>l?3 :4) line  (2),  ray  (3),  or  segment  (4)  intersections 
}//  rvrv/vfv  YAGENAUT@GMAIL  .  COM - LAST-UPDATED-16AUG2014 - 


5.2  Celllntersect()  Parameters 

Z  Z  points  to  a  2-index  array  of  values  that  represents  z(x,y) ,  where 

Z[i][j]  =  z{Xi,yj). 

i  i  specifies  i ,  an  index  that  is  associated  with  the  x  axis.  Values  for  i  cannot  be 

less  than  zero  or  greater  than  m-2,  where  m  is  the  first-index  size  of  Z. 

j  j  specifies  j  ,  an  index  that  is  associated  with  the  y  axis.  Values  for  j  cannot  be 

less  than  zero  or  greater  than  n-2  ,  where  n  is  the  second-index  size  of  Z. 

L  L  specifies  2  points  that  define  a  line  (L={ L0  ,L0s  ,L0  _,Lls^,Lls  ,LX„ }). 

sx  sx  is  the  calculated  value  for  sx ,  the  point  of  intersection  between  L  and  the 

surface  in  scaled  coordinates  (see  Eq.  53).  If  the  function’s  return  value  is  -1  or  0, 
then  sx  is  not  calculated. 
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sy  sy  is  the  calculated  value  for  s  ,  the  point  of  intersection  between  L  and  the 

surface  in  scaled  coordinates  (see  Eq.  53).  If  the  function’s  return  value  is  -1  or  0, 
then  sy  is  not  calculated. 

z  z  is  the  calculated  value  for  z  ,  the  point  of  intersection  between  L  and  the  surface 

(see  Eq.  53).  If  the  function’s  return  value  is  -1  or  0,  then  z  is  not  calculated. 

epsilon  epsilon  specifies  the  minimum  nonzero  value  for  both  a  from  Eq.  50  and  b  from 
Eq.  5 1 .  Thus,  if  |a|  is  less  than  epsilon,  then  a  is  considered  to  be  zero,  and  if  |b| 

is  less  than  epsilon,  then  b  is  considered  to  be  zero. 

5.3  Celllntersect()  Return  Value 

In  addition  to  (possibly)  calculating  values  for  sx,  sy,  and  z,  the  Celllntersect()  function  returns  1 
of  6  values  (-1,  0,  1,2,  3,  or  4): 

•  A  return  value  of-1  indicates  that  a  ,  b  ,  and  c  from  Eq.  28  are  all  effectively  zero.  This 
means  that  line  L  intersects  the  surface  at  an  infinite  number  of  points.  For  a  return  value 
of-1,  sx,  sy,  and  z  are  not  calculated. 

•  A  return  value  of  0  indicates  that  the  line  defined  by  L  does  not  intersect  the  surface  at  any 
point  (and,  thus,  sx,  sy,  and  z  are  not  calculated). 

•  A  return  value  of  1  indicates  that  the  line  defined  by  L  intersects  the  surface  at  a  point  that 
is  outside  the  boundaries  of  the  cell. 

•  A  return  value  of  2  (or  greater)  indicates  that  the  line  defined  by  L  intersects  the  surface 
within  the  bounds  of  the  cell. 

•  A  return  value  of  3  (or  greater)  indicates  that  the  ray  defined  by  L  (from  L0  to  Lx ) 
intersects  the  surface  within  the  bounds  of  the  cell. 

•  A  return  value  of  4  indicates  that  the  line  segment  defined  by  L  intersects  the  surface 
within  the  bounds  of  the  cell. 

5.4  Celllntersect()  Simple  Example 

The  following  example  begins  by  defining  the  surface  that  is  shown  in  Fig.  5.  The  Celllntersect() 
function  is  then  used  to  determine  the  point  of  intersection  between  the  surface  and  a  vertical  line 
located  at  x  =  1 .5  ,  y  =  0.5 . 
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#include  <cstdio>// . printf() 

#include  "ybilinear .  h"// . yBilinear 

int  main(){//<===============A  SIMPLE  EXAMPLE  USING  THE  Celllntersect ( )  FUNCTION 

const  int  m=5,n=3; 

double  Z[m][n]={l,l,0  ,  1,0,0  ,  0,1,0  ,  0,0,0  ,  1,0,1}; 
int  i=yBilinear: :SafeIndex(m,l. 5), j=yBilinear: :SafeIndex(n, .5); 
double  L[6]={1.5, .5, 1,1. 5, .5,0}; 
double  x,y,z; 

int  intersect=yBilinear : : CellIntersect(Z, i, j, L,x,y, z) ; 

if (intersect==-l)printf ( "Intersection  cannot  be  determined . \n" ) ; 

else  if (intersect==0)printf ("The  line  does  not  intersect  the  surface . \n" ) ; 

else  printf( "Intersection  at  x=%.lf,  y=%.lf,  and  z=%.lf .\n",x,y,z); 

^  J  j /v»v/min;iv<\/Y AG  E  N AUT @GMA XL*  C LAST r^j LJ P DAT  E  D 1 6 AUG  2 0  *1 


OUTPUT: 

Intersection  at  x=1.5,  y=0.5,  and  z=0.5. 


5.5  Celllntersect()  Image  Example 

The  following  example  begins  by  defining  the  surface  that  is  shown  in  Fig.  5.  Note  that  the 
surface  is  a  composite  of  8  simple  surfaces  that  are  each  described  by  Eq.  2 1 .  Next,  the 
Celllntersect()  function  is  used  to  calculate  the  point  of  intersection  between  the  simple  surface 
that  corresponds  with  the  i  =  0 ,  j  =  0  cell  and  107  randomly  chosen  lines.  Functions  from  the 
yBmp  namespace,  along  with  the  Rainbow()  function  (presented  in  Section  4.5),  are  used  to 
create  a  pseudo-color  image  of  the  intersections.  Finally,  the  image  is  written  to  a  BMP  file, 
which  is  displayed  in  Fig.  8.  The  color-to-value  key  is  shown  in  Fig.  7.  Since  the  intersecting 
lines  are  randomly  generated,  not  all  points  are  represented.  Non-plotted  points  are  shown  in 
white. 


#include  <cstdlib>// . rand( ) , RANDJMAX 

#include  "y  bmp.h"// . yBmp 

#include  "y  bilinear .  h"// . yBilinear,<cmath>{fabs( )} 

int  main( ){//<==============CREATE  AN  IMAGE  FROM  A  SURFACE  USING  CelllntersectQ 

const  int  m=5, n=3,M=1000,N=500; 

double  Z[m] [n]={l,l,0  ,  1,0,0  ,  0,1,0  ,  0,0,0  ,  1,0,1}; 
unsigned  char*I=yBmp: :NewImage(M,N,255); 
for(int  i=0; i<10000000;++i){ 

double  L[6]={rand( )*(m-l)*l./RAND_MAX, 
rand()*(n-l)*l./RAND_MAX, 

- 5+rand ( )  *10 . /RAND_MAX, 

rand()*(m-l)*l./RAND_MAX, 

rand()*(n-l)*l./RAND_MAX, 

- 5+rand ( )  *10 . /RAND_MAX}; 

double  x,y,z; 

if (yBilinear : : Celllntersect (Z, 0,0, L,x,y,z)>0){ 
int  q=int(x/(m-l)*M) , p=int(y/(n-l)*N); 

if (q>=0&&q<M&&p>=0&&p<N)Rainbow(yBmp : : Get Pixel ( I, q, p) , z,0, 1) ;  }} 
yBmp : :WriteBmpFile( "cell_intersect . bmp",  I); 

^  J  j <v<vi/vivfN»vY AG  E  N AUT (o)GMA XL*  C LAST r^j U P DAT  E  D ^  1 6 AUG  2 0  *1 
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Fig.  8  Image  generated  by  example  code  from  Section  5.5. 


6.  Calculating  Line-Surface  Intersections:  The  Line()  Function 


The  Line()  function  is  designed  to  work  with  the  Celllntersect()  function  to  determine  the  point 
of  intersection  (if  it  exists)  between  a  line  segment,  L  ,  and  an  evenly  spaced  rectangular  surface 
grid.  This  is  accomplished  by  first  using  the  Line()  function  to  find  a  set  of  cells  that  L  may 
intersect,  then  using  the  Celllntersect()  function  to  test  each  of  the  cells  for  intersection. 

The  set  of  cells  specified  by  the  Line()  function  is  typically  much  smaller  than  the  set  of  all  cells 
that  makes  up  a  surface.  Thus,  the  primary  purpose  of  the  Line()  function  is  to  reduce  the  number 
of  cells  that  must  be  tested  using  the  Celllntersect()  function.  For  some  scenarios,  it  may  be 
possible  to  further  reduce  the  number  of  cells  that  need  to  be  tested  by  performing  additional 
tests  after  calling  the  Line()  function  but  before  calling  the  Celllntersect()  function. 

The  Line()  function  uses  the  algorithm  presented  in  Section  2.4  to  determine  the  indices  of  the 
cells  that  line  segment  L  might  intersect  when  traveling  from  L0  to  L] .  However,  the  algorithm 
has  been  optimized  in  2  ways. 

Eq.  36,  which  is  used  to  find  i, ,  can  be  written  to  allow  for  some  calculations  to  be  removed 
from  the  outer  loop: 

xt  =  O’  +  bj  )M  +  P  ,  (54) 


where 
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b 


1  for  A j  =  1 
0  otherwise 


(55) 


for  j  +  J 


otherwise 


(56) 


and 


(57) 


The  Line()  function  has  been  written  to  take  advantage  of  the  fact  that  the  algorithm  presented  in 
Section  2.4  is  more  efficient  for  lines  whose  slopes  are  closer  to  horizontal  than  vertical.  If 
| /  —  /|  >  \  j  —  j\  then  the  algorithm  is  used  as  described.  Otherwise,  the  location  where  the 
projection  of  L  onto  the  x—y  plane  crosses  a  constant*  reference  line  (rather  than  a  constant  y 

reference  line)  is  found.  This  results  in  modified  outer  and  inner  loops  that  are  very  similar  to 
those  in  Section  2.4,  but  where  indices  and  measured  values  associated  with  the  *  and  y  axes 

have  been  exchanged. 

6.1  Line()  Code 

inline  int  Line(//<=====================================LINE-GRID  INTERSECTION 

int  m, int  n,//< . NUMBERS  OF  GRID  INDICES  (m  AND  n  SHOULD  BE  AT  LEAST  2) 

const  double  L[6],//<--A  SCALED  LINE  SEGMENT  {L0SX, L0SY, L0Z, L1SX, L1SY, LIZ} 

int*Aj int*B){//< . STORAGE  FOR  CELL  INDICES  (FOR  EACH,  SIZE  >=  m+n-3) 

int  i=SafeIndex(m, *L) , j=SafeIndex(n, L [ 1 ] ) , I=SafeIndex(m, L [ 3 ] ) , 

D=SafeIndex(n, L[4] ) ,di=i<I?l : -l,d j=j<3  ?1 : -1, bi=di>0?l : 0, bj=d j  >0?1 : 0, k=0; 
double  M=J-j?(L[3]-*L)/(L[4]-L[l]):0,N=I-i?l/M:0,P=*L-L[l]*M,Q=-P*N; 
if (abs(I-i)>abs(3- j))for(int  it;dj*(3- j)>=0; j+=dj,i-=di)//dj*  selects  <  or  > 
for(it=j-D?SafeIndex(m, ( j+bj)*M+P) :I;di*(i-it)<=0; i+=di)A[k]=i,  B[k++]=j; 

else  for(int  jt;di*(I-i)>=0;i+=di, j - =d j ) // . di*  selects  <  or  > 

for( jt=i-I?SafeIndex(nj (i+bi)*N+Q) : 3;dj*( j- jt)<=0; j+=d j)A[k]=i, B[k++]=j; 

return  k;// . the  number  of  elements  to  check  in  A  and  B 

}//  /v/vrv/v  YAGENAUT@GMAIL  .  COM - LAST-UPDATED-16AUG2014 -  | 


6.2  Line()  Parameters 

m  m  specifies  m  ,  the  number  of  horizontal  values  that  the  surface  grid  contains,  m 


L 


n 


must  be  greater  than  1 . 

n  specifies  n  ,  the  number  of  vertical  values  that  the  surface  grid  contains,  n  must 
be  greater  than  1 . 

L  specifies  2  points  that  define  a  line  segment  (L={  L  ,  Lx  _ }). 


B 


A 


A  points  to  storage  for  x-axis  indices.  The  size  of  A  must  be  at  least  m  +  n-3. 
B  points  to  storage  fory-axis  indices.  The  size  of  B  must  be  at  least  m  +  n-3. 
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6.3  Line()  Return  Value 

In  addition  to  calculating  values  for  A  and  B,  the  Line()  function  returns  the  number  of  index 
pairs  that  are  generated.  The  indices  found  in  A  and  B  are  stored  in  the  order  in  which  they  are 
found  when  traveling  from  L0  to  Lx . 

6.4  Line()  Simple  Example 

The  following  example  begins  by  defining  the  surface  that  is  shown  in  Fig.  5.  The  Line()  and 
Celllntersect()  functions  are  then  used  to  determine  the  point  of  intersection  between  the  surface 
and  a  vertical  line  located  at  *  =  1.5 ,  J-0.5. 


#include  <cstdio>// . printfQ 

#include  "ybilinear .  h"// . yBilinear 

int  main(){//<========================A  SIMPLE  EXAMPLE  USING  THE  LineQ  FUNCTION 

const  int  m=5,n=3; 

double  x0=0,y0=0,dx=l,dy=l,Z[m] [n]={l,l,0  ,  1,0,0  ,  0,1,0  ,  0,0,0  ,  1,0,1}; 
double  L[6]={1.5, .5, 1,1. 5, .5,0}; 

int  *A=new  int[m+n] , *B=new  int [m+n] , k=yBilinear : : Line(m, n ,  L, A, B) ; 
double  x,y,z; 
bool  b=0; 

for(int  j=0,q; j<k;++j) 

if ( (q=yBilinear : :CellIntersect(Z, A[ j ] ,B[j],L,x,y,z))>2){ 
b= ! (q-4) ; break; } 

printf(  "Intersection  at  x=%.lf,  y=%.lf,  and  z=%.  If .  \n",x,y,  z); 

J  j  /  (m/v/n/iv/n/ivY  AG  E  N  AUT  (o)GMA  XL*  C  LAST  r^j  LJ  P  DAT  E  D  r^j  1 6  AUG  2  0 1 


OUTPUT: _ 

Intersection  at  x=1.5,  y=0.5j  and  z=0.5. 


6.5  Line()  Image  Example 

The  following  example  begins  by  defining  the  surface  that  is  shown  in  Fig.  5.  Next,  the  Line() 
and  Celllntersect()  functions  are  used  to  calculate  the  point  of  intersection  between  the  surface 
and  107  randomly  chosen  line  segments.  Functions  from  the  yBmp  namespace,  along  with  the 
Rainbow()  function  (presented  in  Section  4.5),  are  used  to  create  a  pseudo-color  image  of  the 
intersections.  Finally,  the  image  is  written  to  a  BMP  file,  which  is  displayed  in  Fig.  9.  The  color- 
to-value  key  is  shown  in  Fig.  7.  Since  intersecting  lines  are  randomly  generated,  not  all  points 
are  represented.  Non-plotted  points  are  shown  in  white. 


#include  <cstdlib>// . rand( ) ,  RAND  MAX 

#include  "y  bmp.h"// . yBmp 

#include  "y  bilinear .  h"// . yBilinear,<cmath>{fabs( )} 

int  main( ){//<=======================CREATE  AN  IMAGE  FROM  A  SURFACE  USING  LineQ 

const  int  m=5, n=3JM=1000JN=500; 

double  Z[m][n]={l,l,0  ,  1,0,0  ,  0,1,0  ,  0,0,0  ,  1,0,1}; 
unsigned  char*I=yBmp: :NewImage(M,N,255); 
int  K=0, 3=0, *A=new  int[m+n],*B=new  int[m+n]; 
for(int  i=0; i<10000000;++i){ 

double  L[6]={rand( )*(m-l)*l./RAND_MAX, 
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rand()*(n-l)*l./RAND_MAX, 

- 5+rand ( )  *10 . /RAND_MAX, 

rand()*(m-l)*l./RAND_MAX, 

rand()*(n-l)*l./RAND_MAX, 

- 5+rand ( )  *10 . /RAND_MAX}; 

double  XjYjZ; 

int  k=yBilinear : : Line(m, n, L,Aj B) ; 
bool  b=0; 

for(int  j=0,q; j<k;++j) 

±-F( (q=yB±linear : : CellIntersect(Z,A[ j ]JB[j]JLJxJy,z))>2){ 
b= ! (q-4) ; break; } 
if(b){ 

int  j=int(x/(m-l)*M) , k=int(y/(n-l)*N); 

if ( j>=0&&j<M&&k>=0&&k<N)Rainbow(yBmp : : Get Pixel (I , j, k) , z,0, 1) ; }} 
yBmp : :WriteBmpFile( "intersect . bmp" , I) ; 

^  j  j /v<v»ni/v/\hv Y /\(j  E  N AUT(g)GMAI  L  •  C L /\ S T ^ LJ P D/\”F  E  D ^  1 6 /\ U G 201 r^j 


Fig.  9  Image  generated  by  example  code  from  Section  6.5. 

Note  that  the  bottom-left  comer  of  the  image  is  similar  to  the  image  that  is  displayed  in  Fig.  8. 
However,  the  image  from  Fig.  8  has  fewer  missing  pixels  than  the  image  from  Fig.  9.  This  is 
because  different  intersection  criteria  were  used  to  generate  the  2  images;  only  line-segment 
intersections  were  plotted  in  Fig.  9. 
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7.  Calculating  Surface  Gradients:  The  Gradient^)  Function 


The  Gradient()  function  uses  Eq.  39  to  calculate  the  surface  gradient  at  a  user-specified  location. 
To  simplify  the  Gradient()  function,  Eq.  39  has  been  rewritten  in  terms  of  scaled  coordinates 
using  the  substitutions  from  Eq.  44: 


Eq.  60  can  be  used  to  find  (V  </>)x  and  (V  <j>)  r : 

cm 

Ax 


ZMj)~(Sy  -J)(Zi,J+l 

(58) 

ziJ+i)-(sx-i)(zMJ  -zMJ+1) . 

(59) 

cm 

and  (V</>)v  - 

(60) 

Ay 

The  components  of  the  gradient  can  be  used  to  determine  the  direction  of  steepest  ascent,  the 
slope  in  the  direction  of  steepest  ascent  (Eq.  40),  and  a  unit  vector  that  is  normal  to  the  surface 
(Eq.  41). 

7.1  Gradient()  Code 

template<class  T>void  Gradient(//< . SURFACE  GRADIENT 

const  T&Z,//< . -POINTER  TO  Z  WHERE  Z[i]  [  j]=Z(X[i] ,Y[  j] ) 

int  i,int  j,//< . GRID  INDICES  (0<=i<m-l  &  0<=j<n-l) 

double  sx, double  sy ,//< . A  GRID  LOCATION  IN  SCALED  COORDINATES 

double&delsx,double&delsy){//< . SCALED  GRADIENT  COMPONENTS  (CALCULATED) 

delsx=(sy- j-l)*(Z[i] [ j ] -Z[i+1] [j])-(sy-j)*(Z[i] [ j+1] -Z[i+1] [ j+1] ) J 
delsy=(sx-i-l)*(Z[i] [j] -Z[i] [j+1] ) - (sx-i)*(Z[i+l] [ j ] -Z[i+1] [ j+1] ) ; 

}//  Y AG  E  N AUT (3GMAI L  .  L AS T ~UPDAT  ED~16AUG2014-~~~~~~ 


7.2  Gradient^)  Parameters 

Z  Z  points  to  a  2-index  array  of  values  that  represents  z(x,y) ,  where 

Z[i][j]  =  z(xi,yj). 

i  i  specifies  i,  an  index  that  is  associated  with  the  x  axis.  Values  for  i  cannot  be 

less  than  zero  or  greater  than  m-2,  where  m  is  the  first-index  size  of  Z. 

j  j  specifies  j  ,  an  index  that  is  associated  with  the  y  axis.  Values  for  j  cannot  be 

less  than  zero  or  greater  than  n  -  2  ,  where  n  is  the  second-index  size  of  Z. 

sx  sx  specifies  sx  ,  a  scaled  distance  from  Eq.  42. 

sy  sy  specifies  s  ,  a  scaled  distance  from  Eq.  42. 
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delsx  delsx  is  the  calculated  value  for  (V  <f>)  s  ,  the  x-component  of  the  surface  gradient 

(in  scaled  coordinates)  from  Eq.  58. 

delsy  delsy  is  the  calculated  value  for  (V^)s  ,  the  v-component  of  the  surface  gradient 

(in  scaled  coordinates)  from  Eq.  59. 

7.3  Gradient()  Simple  Example 

The  following  example  begins  by  defining  the  surface  that  is  shown  in  Fig.  5.  The  Gradient() 
function  is  then  used  to  estimate  the  slope  of  the  surface  at  x  =  1 .5  and  y  =  0.5 . 


#include  <cstdio>// . printfQ 

#include  "y  bilinear . h"// . yBilinear, <cmath>{sqrt( )} 

int  main(){//<====================A  SIMPLE  EXAMPLE  USING  THE  Gradient()  FUNCTION 

const  int  m=5,n=3; 

double  Z[m][n]={l,l,0  ,  1,0,0  ,  0,1,0  ,  0,0,0  ,  1,0,1}; 
double  x=1.5,y=.5; 

int  i=yBilinear: :SafeIndex(m,x), j=yBilinear: :SafeIndex(n,y); 
double  delx,dely;/*<-*/yBilinear : : Gradient (Z, i, j,x,y,delx,dely); 
printf("At  x=%.lf  and  y=%.lf,  slope=%.lf\n",x,y,sqrt(delx*delx+dely*dely)); 

^  I  /  AG  E  N  AUT  (o)GMA  XL*  C  LAST  r^j  LJ  P  DAT  E  D rKj  1 6  AUG  2  0  X  4^  ^  ~  ~  ~  ^ 


OUTPUT: 

At  x=1.5  and  y=0.5,  slope=0.0 


7.4  Gradient()  Image  Example 

The  following  example  begins  by  defining  the  surface  that  is  shown  in  Fig.  5.  Next,  functions 
from  the  yBmp  namespace  are  used,  along  with  the  Gradient()  function  and  the  Rainbow() 
function  (presented  in  Section  4.5),  to  create  a  pseudo-color  image  of  the  gradient  of  the  surface. 
Finally,  the  image  is  written  to  a  BMP  file,  which  is  displayed  in  Fig.  10.  The  color-to-value  key 
is  shown  in  Fig.  7. 

Note  that  at  the  boundaries  between  cells,  the  gradient  may  contain  discontinuities.  The 
discontinuities  are  the  result  of  the  manner  in  which  the  surface  has  been  constructed. 


#include  "y  bmp.h"// . yBmp 

#include  "y  bilinear .  h"// . yBilinear,<cmath>{fabs( ) ,  sqrt( )} 

int  main( ){//<===================CREATE  AN  IMAGE  FROM  A  SURFACE  USING  GradientQ 

const  int  m=5, n=3,M=1000,N=500; 

double  Z[m] [n]={l,l,0  ,  1,0,0  ,  0,1,0  ,  0,0,0  ,  1,0,1}; 
unsigned  char*I=yBmp: :NewImage(M,N,255); 
for(int  p=0;p<M;++p)for(int  q=0;q<N;++q){ 

double  x=p*(m-l)*l./(M-l) ,y=q*(n-l)*l./(N-l); 

int  i=yBilinear: :SafeIndex(m,x), j=yBilinear: :SafeIndex(n,y); 

double  delx,dely;/*<-*/yBilinear : : Gradient (Z, i, j,x,y,delx,dely); 

Ra inbow (yBmp : :GetPixel(I, p,q) , sqrt(delx*delx+dely*dely) ,0, 1.5);} 
yBmp : :WriteBmpFile ("gradient. bmp", I); 

"^11  ‘'"'"'"'"'"''YAGENAUT^IGMAI  L  •  C 'W rv rw rv rv rv rw rv rw rv rw rv rw /V) rw rv rw rv  L /\ ^ T ^ U  P  D/\T  E  D ^  1 0 /\ U G  201 4^^ ^ ^ 
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Fig.  10  Image  generated  by  example  code  from  Section  7.4. 


8.  Performance 


The  following  example  estimates  the  calculation  time  per  107  iterations  for  the  Interpolate!), 
Celllntersect(),  and  Gradient!)  functions  by  randomly  sampling  and  intersecting  a  randomly 
generated  surface.  The  yRandom  namespace  is  used  to  generate  pseudorandom  numbers.  Figure 
1 1  summarizes  the  results.  Time  values  will  vary  based  on  computer  specifications,  compiler, 
compiler  settings,  etc. 


#include  <cstdio>// . printfQ 

#include  <ctime>// . clock( ) ,  CLOCKS_PER_SEC 

#include  "ybilinear .  h"// . yBilinear 

#include  "y_random. h"// . yRandom 

int  main( ){//<=TEST  THE  SPEEDS  OF  Interpolate! ) j  CelllntersectQ,  and  Gradient!) 

const  int  n=14jN=l<<njM=10000000; 

unsigned  I [625] ; /*<-*/yRandom: : Initialize^, 1) state  of  Mersenne  twister 
double**Z=new  double*[N];/*<-*/for(int  i=0; i<N;++i)Z[i]=new  double[N]; 
for (int  i=0; i<N;++i)for(int  j=0; j<N;++j)Z[i] [j]=y Random: :RandU(I,0,l); 
printf("  size  |  Interpolation!)  I  Celllntersect( )  |  Gradient!" 

")\n  of  | . | . -| . " 

" . \n  array  |  |  z  |  |  intersect  |  |  " 

"  delx\n  (m)  |  time  (s)  |  avg.  |  time  (s)  |  avg.  |  time  (s)  |" 

”  avg. \n  . -  | . -  | . | . | . | . " 

"| . \n")j// . table  header 

for(int  m=2;m<=N;m*=2){ 

double  s=0jt=clock( ) ; // . begin  test  for  Interpolate!) 

for(int  k=0; k<M;++k){ 

double  x=yRandom: :Randll(I,0,m-l),y=yRandom: :RandU(I,0,m-l); 
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int  i=yBilinear: :SafeIndex(m,x), j=yBilinear: :SafeIndex(m,y); 

double  z=yBilinear: :Interpolate(ZJiJ j,x,y); 

s+=z;} 

printf ( "%7d  |%8.3f  |%8.3f  | ",m, (clock( ) -t)/CLOCKS_PER_SECJ s/M) ; 

s=0jt=clock( ) ; // . begin  test  for  Celllntersect( ) 

for(int  k=0; k<M;++k){ 

double  L[6]={yRandom: :RandU(IJ0Jm-l) ,yRandom: :RandU(IJ0^m-l)J 
y Random: :  Randll(I,  -1,  l),y Random: :  Randll(1 , 0jm-l) , 
y Random: : RandU(I j0,m-l) ,y Random: : RandU(I , -1,1) }; 
double  x,y,z; 

int  i=yRandom: :RandI(IJ0Jm-2)J j=yRandom: :RandI(IJ0Jm-2); 
int  test=yBilinear : : CellIntersect(Z, i, j , L,x,y,z) ; 
if (test==4)s++; } 

printf ( "%9 . 3f  |%8.3f  | ”, (clock( ) -t)/CLOCKS_PER_SECJ s/M) ; 

s=0jt=clock( ) ; // . begin  test  for  GradientQ 

for(int  k=0; k<M;++k){ 

double  x=yRandom: :RandU(Ij0,m-l) ,y=yRandom: :RandU(IJ0Jm-l)JdelxJdely; 
int  i=yBilinear: :SafeIndex(m,x), j=yBilinear: : Safelndex(m,y) ; 
yBilinear : :Gradient(Z, i, j,Xjyjdelx,dely); 
s+=delx; } 

printf ("%9.3f  |%8 . 3f \n", (clock( ) -t)/CLOCKS_PER_SECJ s/M) ; } 

^  J  j AG  E  N AUT (3GF1A XL*  C LAST ^ LJ P DAT  E  D ^  1 6 AUG  2 0 1 


OUTPUT: 


size  | 
of  | 

array  j 
(m)  j 

Interpolation( )  | 

i 

Celllntersect( )  | 

i 

Gradient( ) 

1 

time  (s)  | 
_ 1 

1 

Z  1 

avg. 

_ | 

1 

time  (s)  | 
_ 1 

i 

intersect | 
avg.  | 

_ 1 

1 

time  (s)  | 
_ 1 

delx 

avg. 

2  i 

1 

0.234  | 

1 

0.672  | 

1 

1.201  | 

1 

0.275  | 

1 

0.234  | 

0.071 

4  | 

0.297  | 

0.439  j 

1.216  j 

0.055  j 

0.297  j 

0.201 

8  | 

0.296  j 

0.439  j 

1.248  j 

0.014  j 

0.297  j 

-0.020 

16  j 

0.280  | 

0.486  j 

1.233  j 

0.005  | 

0.265  j 

-0.008 

32  | 

0.281  j 

0.494  | 

1.232  j 

0.002  | 

0.281  j 

-0.001 

64  | 

0.281  | 

0.505  | 

1.263  j 

0.001  | 

0.281  j 

-0 . 000 

128  j 

0.297  j 

0.504  | 

1.357  j 

0.001  j 

0.296  j 

-0.000 

256  j 

0.312  i 

0.501  | 

1.420  j 

0 . 000  | 

0.296  j 

0 . 000 

512  | 

0.328  i 

0.500  j 

1.404  j 

0 . 000  j 

0.296  j 

-0.000 

1024  | 

0.484  | 

0.500  | 

1.591  j 

0 . 000  | 

0.452  j 

-0.000 

2048  | 

0.905  j 

0.500  j 

2.028  j 

0 . 000  | 

0.874  j 

0 . 000 

4096  j 

1.014  j 

0.500  j 

2.106  j 

0 . 000  | 

1.045  j 

0 . 000 

8192  j 

1.154  j 

0.500  j 

2.372  j 

0 . 000  j 

1.076  j 

-0.000 

16384  j 

1.451  j 

0.500  j 

2.558  j 

0 . 000  j 

1.311  j 

-0.000 
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InterpolationQ  Performance 
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Fig.  11  Average  time  per  iteration  for  the  InterpolateQ,  CelllntersectQ,  and  GradientQ  functions 


9.  Performance  II 


The  following  example  begins  by  defining  2  functions  that  can  be  used  to  find  intersections 
between  lines  and  evenly  spaced  rectangular  surface  grids.  The  first  function,  Simplelntersect(), 
uses  the  Celllntersect()  function  to  check  all  of  the  grid’s  cells.  The  second  function, 
Efficientlntersect(),  uses  the  Celllntersect()  function  to  check  a  subset  of  the  grid’s  cells,  as 
specified  by  the  Line()  function. 

Next,  the  example  code  estimates  the  calculation  time  per  105  iterations  for  the  Simplelntersect() 
and  Efficientlntersect()  functions.  The  yRandom  namespace  is  used  to  generate  pseudorandom 
numbers. 
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The  example  demonstrates  2  things:  First,  that  using  the  Line()  function  to  reduce  the  number  of 
cells  to  check  for  intersection  can  result  in  large  increases  in  performance.  Second,  average 
intersect  values  match  between  the  2  methods,  providing  evidence  that  the  algorithm  for  the 
Line()  function  is  valid. 


template<class  T>bool  Simplelntersect(//<=================LINE-GRID  INTERSECTION 

const  T&Z,  //< -  - . - . -POINTER  TO  Z  WHERE  Z[i]  [  j  ]=Z(X[i]  ,Y[  j  ]  ) 

int  m, int  n,//< . GRID  SIZES 

const  double  L  [6] , //< . A  LINE  SEGMENT  {  L0X,  L0Y,  L0Z,  L1X,  L1Y,  LIZ} 

double&sx,double&sy,double&sz, //< . --POINT  OF  INTERSECTION  (CALCULATED) 

double  epsilon=lE-9){//< . -VALUE  FOR  DIVISION-BY-ZERO  CHECK 

double  d=-lJdtJxt,ytJzt; 

for(int  i=0; i<m-l;++i)for(int  j=0; j<n-l;++j) 

if (yBilinear : :CellIntersect(Z, i, j, L,xt ,yt , zt , epsilon )==4){ 
xt-=L[0],yt-=L[l],zt-=L[2]; 
if (dt=xt*xt+yt*yt+zt*zt<d | |d<0) 

d=dt, sx=xt+L[0] , sy=yt+L[l] , sz=zt+L[2] ; } 
return  d<0?0:l; 

^  j  j AG  E  N AUT (o)GMA XL*  C LAST ^ LJ P DAT  E  D ^  1 6 AUG  2 0  X 4^^ ^ ^ 


templatecclass  T>bool  Eff icientIntersect(//<==============LINE-GRID  INTERSECTION 

const  T&Zj  //< . POINTER  TO  Z  WHERE  Z[ i]  [  j  ]  =Z(X[i] ,  Y[  j  ] ) 

int  m, int  n,//< . GRID  SIZES 

const  double  L[6],//< . A  LINE  SEGMENT  { L0X,  L0Y,  L0Z,  L1X,  L1Y,  LIZ} 

double&sx,double&sy,double&sz,//< . POINT  OF  INTERSECTION  (CALCULATED) 

int*A,int*B,//< . STORAGE  FOR  CELL  INDICES  (FOR  EACH,  SIZE  >=  m+n-3) 

int  k,//< . RETURN  VALUE  FROM  THE  Line()  FUNCTION 

double  epsilon=lE-9){//< . VALUE  FOR  DIVISION-BY-ZERO  CHECK 

int  q=0;/*<-*/for(int  j=0; j<k;++j) 

if ( (q=yBilinear : :CellIntersect(Z,A[ j ] , B[ j ] , L, sx, sy, sz, epsilon ) ) >2) break; 
return  !(q-4); 

"^11  ‘'"'"'"'"'"''YAGENAUT^IGMAI  L  •  C  01^^^  'v,v  ^  ~  ^  ~  ~  ^  rv  rv  rv  rv  rw  rv  rv  rv  rsj  rv  rw  rv  L  A  S  T  rv/  U  P  D/\T  E  D  ^  1 0  /\  U  G  201 4^^  ^ 


#include  <cstdio>// . printf() 

#include  <ctime>// . clock( ) ,CLOCKS_PER_SEC 

#include  "ybilinear .  h"// . yBilinear 

#include  "y_random.  h"// . y Random 

int  main( ){//<========C0MPARE  LineQ  FUNCTION  PERFORMANCE  TO  AN  EXAUSTIVE  METHOD 

const  int  n=ll,N=l<<n,M=100000; 

unsigned  I [625] ;/*< -*/yRandom: : Initialize^, 1) state  of  Mersenne  twister 
double**Z=new  double*[N];/*<-*/for(int  i=0; i<N;++i)Z[i]=new  double[N]; 
for (int  i=0; i<N;++i)for(int  j=0; j<N;++j)Z[i] [j]=y Random: :RandU(I,0,l); 


printfC 


size 

of 

array 

(m) 


Simplelntersect( ) 


time  (s) 


intersect 

avg. 


Ef f icient Inter sect ( )\n" 

. \n’ 

|  intersect\n" 
time  (s)  |  avg.\n" 


An"); 


for(int  m=2;m<=N;m*=2){ 
double  s,t; 

s=0,t=clock( ) ,yRandom: : Initialize (1 , 1) ; 
for(int  k=0; k<M;++k){ 
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double  L[6]={yRandom: :RandU(IJ0Jm-l)JyRandom: :Randll(I,0,m-l), 
y Random: : Randll(I, -1, 1) ,y Random: : RandU(I , 0,m-l) , 
y Random: : RandU(I ,0,m-l) ,y Random: : Randll(I, -1,1) }; 
double  x,y,z; 

int  i=yRandom: :RandI(I,0,m-2), j=yRandom: :RandI(I,0,m-2); 
bool  test=SimpleIntersect(Z,m,m, L,x,y,z); 
if (test)s++; } 

printf ( "%7d  |%9.3f  |%10.7f  | ",m, (clock( ) -t)/CLOCKS_PER_SEC, s/M) ; 

s=0,t=clock( ) ,yRandom: : Initialize^,  1) ; 
int  *A=new  int[m+m-3] , *B=new  int[m+m-3]; 
for(int  k=0; k<M;++k){ 

double  L[6]={yRandom: :RandU(I,0,m-l),yRandom: :RandU(I,0,m-l), 
y Random: : RandU(I, -1, l),y Random: : RandU(I , 0,m-l) , 
y Random: : RandU(I ,0,m-l) ,y Random: : RandU(I, -1,1)}; 
int  K=yBilinear: : Line(m,m, L, A, B) ; 
double  x,y,z; 

int  i=yRandom: :RandI(I,0,m-2), j=yRandom: :RandI(I,0,m-2); 
bool  test=Efficient Intersect (Z,m,m, L,x,y,z,A, B, K); 
if (test)s++; } 
delete [] A, delete [ ]B; 

printf ("%9.3f  |%10.7f \n" , (clock( ) -t)/CLOCKS_PER_SEC, s/M) ; } 

^  I  j  AG  E  N  AUT  (3GF1A  XL*  C  LAST  r^j  LJ  P  DAT  E  D  ^  1 6  AUG  2  0 1 


OUTPUT: 


size 

|  Simplelntersect( ) 

1 

Eff icient Inter sect ( ) 

of 

1 _ 

1 

1 

array 

1 

intersect 

1 

1 

intersect 

(m) 

|  time  (s) 
_i _ 

1 

I 

avg. 

1 

1 

time  (s) 

1 

I 

avg. 

2 

1 

|  0.015 

1 

1 

0.2961400 

1 

1 

0.015 

1 

1 

0.2961400 

4 

j  0.047 

1 

0.2324000 

1 

0.016 

1 

0.2324000 

8 

j  0.219 

1 

0.4265600 

1 

0.047 

1 

0.4265600 

16 

j  1.061 

1 

0.5226000 

1 

0.046 

1 

0.5226000 

32 

|  4.166 

1 

0.5722800 

1 

0.078 

1 

0.5722800 

64 

|  17 . 044 

1 

0.60912 00 

1 

0.110 

1 

0.6091200 

128 

j  67.985 

1 

0.6480700 

1 

0.218 

1 

0.6480700 

256 

j  272.998 

1 

0.6754100 

1 

0.406 

1 

0.6754100 

512 

j  1089.319 

1 

0.6941000 

1 

0.748 

1 

0.6941000 

1024 

j  4402.322 

1 

0.7071800 

1 

1.498 

1 

0.7071800 

2048 

j 17427.493 

1 

0.7174700 

1 

3.822 

1 

0.7174700 

10.  Surface-Elevation  Example:  The  Far  Side  of  the  Moon 


The  image  shown  in  Fig.  12  was  obtained  from  NASA’s  Website.  It  presents  a  topographic  view 
of  the  far  side  of  the  moon.  Prior  to  working  with  the  image,  the  file  was  converted  from  JPG 
format  to  BMP  format  using  Microsoft  Paint. 
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Reprinted  from  NASA’s  Goddard  Space  Flight  Center/DLR/ASU.  LRO  Camera  Team  Releases  High  Resolution  Global 
Topographic  Map  of  Moon,  16  Nov  2011.  [accessed  2014  May  29].  http://www.nasa.gov/mission_pages/LRO/news 
/lro-topo.html. 

Fig.  12  Topographic  image  of  the  far  side  of  the  moon 

The  example  code  that  follows  begins  by  using  the  yBmp  namespace  to  read  the  image  file  that 
is  presented  in  Fig.  12  into  memory.  Next,  the  code  uses  a  vertical  line  of  pixels  to  interpret  the 
scale  that  is  shown  to  the  right  of  the  moon.  The  vertical  line  is  placed  horizontally  in  the  center 
of  the  scale,  with  the  starting  location  at  the  bottom-most  position  in  the  scale  and  the  ending 
location  at  the  top-most  position. 

Next,  the  code  calculates  elevation  values  for  a  portion,  outlined  in  black,  of  the  image  of  the 
moon  in  Fig.  12.  This  is  accomplished  by  reading  pixel  values  from  the  image  and  performing 
nearest-neighbor  searches  between  the  pixel  color  values  and  a  path  in  3-dimensional  Cartesian 
color  space  that  is  defined  by  the  values  that  have  been  read  from  the  scale.  Although  the  scale  in 
figure  12  isn’t  linear,  for  this  simple  example,  it’s  treated  as  if  it  is. 

Finally,  the  code  creates  3  images  based  on  the  elevation  data  obtained  from  Fig.  12. 


Elevation  (m) 

—  ^  0760 


-8769 


-T2796 


->1186 


- -5168 


9150 


LROC  WAC  Topography  80  S  to  80  N 
LOLA  80°N,S  to  the  poles 


Orthographic  projection  centered 
on  the  farside 
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#include  "y_bmp.h"// . yBmp 

#include  "ybilinear . h"// . yBilinear,<cmath>{fabs( ), sqrt( )} 

int  main( ){//<==============CREATE  VARIOUS  IMAGES  BASED  ON  MOON  TOPOGRAPHIC  DATA 

// . READ  ORIGINAL  IMAGE  AND  CAPTURE  SCALE- . 

unsigned  char*I=yBmp : : ReadBmpFile( "604359main_WAC_CSHADE_O000N1800_1000. bmp" ); 

double*X=new  double[582];// . scale  values  [0,1] 

unsigned  char**C=new  unsigned  char* [582] ;//. scale  colors  ({B,G,R}  and  [0,255]) 
for (int  k=0; k<582;++k)X[k]=k/581 . ,C[k]=yBmp : :GetPixel(I, 1060, k+143) ; 

// . CAPTURE  AND  INTERPRET  A  PORTION  OF  THE  ORIGINAL  IMAGE . 

double  Z[600] [300] ; // . elevation  values 

for(int  p=0;p<600;++p)for(int  q=0;q<300;++q){ 
unsigned  char*c=yBmp: :GetPixel(I , p+200, q+400) ; 
double  d,D=-l; 
for(int  k=0; k<582;++k){ 

d=(C[k][0]-c[0])*(C[k][0]-c[0])+(C[k][l]-c[l])*(C[k][l]-c[l]) 
+(C[k][2]-c[2])*(C[k][2]-c[2]); 
if (d<D | | D<0)D=d,Z[p] [q]=X[k] ; }} 

// . -CREATE  A  NEW  (MONOCHROME)  IMAGE  BASED  ON  ELEVATION  VALUES- . 

unsigned  char*I2=yBmp : :NewImage(600,300,255); 
for(int  p=0;p<600;++p)for(int  q=0;q<300;++q){ 
unsigned  char*c=yBmp: :GetPixel(I2, p, q) ; 
c[0]=c[l]=c[2]=int(Z[p][q]*255);} 
yBmp : :WriteBmpFile( "moon_topography . bmp",  12) ; 

// . CREATE  A  SURFACE -OCCULTATION  IMAGE . 


double  L[6]={0,0,5};// . vantage  point 

int*A=new  int[600+300-3] , *B=new  int[600+300-3] ; 
for(int  p=0;p<600;++p)for(int  q=0;q<300;++q){ 

L[3]=p,  L[4]=q,  L[5]=Z[p]  [q ]  +  . 00000001; // . viewable  point? 

int  K=yBilinear: : Line(600, 300, L, A, B) ,Q=0; 
double  x,y,z; 


for (int  k=0; k<K&&Q! =4;++k)Q=yBilinear : : CellIntersect(Z, A[k] , B[k] , L,x,y, z) ; 
if (Q==4){ 

unsigned  char*c=yBmp: :GetPixel(I2,p,q); 
c[0]=c[l]=c[2]=0; }} 

yBmp : :WriteBmpFile( "moon_occultation . bmp" ,  12) ; 

// . CREATE  SURFACE-GRADIENT  IMAGE- . 

for(int  p=0;p<600;++p)for(int  q=0;q<300;++q){ 

int  i=yBilinear: :SafeIndex(600,p), j=yBilinear: :SafeIndex(300,q); 
double  delx,dely;/*<-*/yBilinear : : Gradient (Z, i, j,p,q,delx,dely); 

Ra inbow (yBmp : :GetPixel(I2, p, q) , sqrt(delx*delx+dely*dely) , 0,  .4);  } 
yBmp : :WriteBmpFile( "moon_gradient . bmp", 12) ; 

^  I  j  AG  E  N  AUT  (3GF1A  XL*  C  LAST  /v,U  P  DAT  E  D  ^  1 6  AUG  2  0 1  ^  ^  r^jr^j 
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Figure  13  presents  a  topographic  view  of  the  moon  with  elevations  denoted  by  shades  of  gray. 


10760  nH 


805  m— 


-9150m 


Fig.  13  Topographic  view  of  a  portion  of  the  far  side  of  the  moon 

Figure  14  was  created  by  placing  a  vantage  point  at  the  lower-left  comer  of  the  image  and  at  an 
elevation  of  90,400  m.  The  Line()  and  Celllntersect()  functions  were  used  to  detennine  line  of 
sight.  Any  point  that  was  not  visible  from  the  vantage  point  was  recolored  to  black. 


Figure  15  shows  the  magnitude  of  the  scaled  surface  gradient  for  the  topographic  information 
that  is  shown  in  Fig.  13. 


Fig.  15  View  of  the  magnitude  of  the  scaled  surface  gradient  for  a  portion  of  the  far  side  of  the  moon 


11.  Code  Summary 


A  summary  sheet  is  provided  at  the  end  of  this  report.  It  presents  the  yBilinear  namespace,  which 
contains  the  SafelndicesQ,  InterpolateQ,  CelllntersectQ,  Line(),  and  GradientQ  functions. 
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y_bilinear. h 


#ifndef  YBILINEARGUARD//  See  Yager,  R.D.  "Working  with  Evenly  Spaced, 

ttdefine  YBILINEARGUARD//  Rectangular  Surface  Grids  Using  C++" 

#include  <cmath>// . fabs(),sqrt() 

#include  <cstdlib>// . abs() 

namespace  yBilinear{// 

inline  int  fg|eIndex(//<=====CALCULATE  A  SAFE  INDEX  AT  A  GIVEN  SCALED  LOCATION 

int  m,//< - NUMBER  OF  GRID  INDICES  (m  SHOULD  BE  AT  LEAST  2) 

double  s){//< - A  GRID  LOCATION  IN  SCALED  COORDINATES 

return  s<0?0:s>m-2?m-2:int(s); 

}// - YAGENAUT@GMAIL .  COM - - LAST~UPDATED~16AUG2014 - 

templatecclass  T>double  lflt#$Qll(t;ei(//<=================BILINEAR  INTERPOLATION 

const  T&Z, / / < - POINTER  TO  Z  WHERE  Z[i]  [ j ]=Z(X[i] ,Y[  j]  ) 

int  i,int  j,//< - GRID  INDICES  (0<=i<m-l  &  0<=j<n-l) 

double  sx, double  sy){//< - A  GRID  LOCATION  IN  SCALED  COORDINATES 

double  za=(sy-j)*(Z[i] [ j+1] -Z[i] [ j] )+Z[i] [ j]; 

return  (sx-i)*((sy-j)*(Z[i+l][j+l]-Z[i+l][j])+Z[i+l][j]-za)+za; 

}// - YAGENAUT@GMAIL .  COM - - LAST~UPDATED~16AUG2014 - 

template<class  T>int  Celllntersect(//<==================LINE-CELL  INTERSECTION 

const  T&Z, / / < - - - POINTER  TO  Z  WHERE  Z[i]  [ j ]=Z(X[i] ,Y[  j]  ) 

int  i,int  j,//< - GRID  INDICES  (0<=i<m-l  &  0<=j<n-l) 

const  double  L[6],//< - A  SCALED  LINE  {L0SX, L0SY, L0Z, L1SX, L1SY, LIZ} 

double&sx, double&sy, double&z, //< - SCALED  INTERSECTION  POINT  (CALCULATED) 

double  epsilon=lE-9){//< - VALUE  FOR  DIVISION-BY-ZERO  CHECK 

double  C0=L[3]-L[0] ,C1=L[4] -L[1],C2=L[5]-L[2]; 
double  C5=i - L [0] , C6= j -L [1] , C3=C5+1,C4=C6+1; 
double  a=C0*Cl*(Z[i] [ j] -Z[i+1] [ j] -Z[i] [ j+l]+Z[i+l] [j+1]), 
b=- (Cl*C3+C0*C4)*Z[i] [ j]+(Cl*C5+C0*C4)*Z[i+l] [ j] 

+(Cl*C3+C0*C6)*Z[i] [ j+1]- (Cl*C5+C0*C6)*Z[i+l] [ j+1] -C2, 
c=C3*C4*Z[i] [ j] -C5*C4*Z[i+l] [ j] -C3*C6*Z[i] [ j+l]+C5*C6*Z[i+l] [j+1] -L[2] ; 
double  t,D; 

if (fabs (a )<epsilon&&fabs(b)<epsilon) return  fabs(c)<epsilon?-l:0; 
if (fabs(a)<epsilon)t=-c/b,sx=C0*t+L[0],sy=Cl*t+L[l],z=C2*t+L[2]; 
else{ 

if((D=b*b-4*a*c)<0)return  0;//..no  intersection,  <sx,sy,sz>  not  calculated 
t=(-b-sqrt(D) )/2/a, sx=C0*t+L[0] , sy=Cl*t+L[l] ,z=C2*t+L[2] ; 
double  tb=(-b+sqrt(D) )/2/a, sxb=C0*tb+L[0],syb=Cl*tb+L[l],zb=C2*tb+L[2]; 
if (t>=0&&t<=l&&tb>=0&&tb<=l){ 

bool  ra=sx<i| |sx>i+l| |sy<j | |sy>j+l?l:0; 
bool  rb=sxb<i| |sxb>i+l| | syb< j | |syb>j+l?l:0; 
if (ra&&rb&&t>tb){sx=sxb,sy=syb,z=zb; return  1; } 
else  if(ra| | !rb&&t>tb)t=tb,sx=sxb,sy=syb,z=zb;} 
else  if (t>tb&&tb>0| |t<=0&&tctb)t=tb,sx=sxb,sy=syb,z=zb;} 

if  (sx<i  |  |  sx>i+l  |  |  sy<j  |  |  sy>j+l)return  1;// . out  of  bounds 

return  t<0?2: (t>l?3:4);//. . . .line  (2),  ray  (3),  or  segment  (4)  intersections 

}// - YAGENAUT@GMAIL .  COM - LAST-UPDATED-16AUG2014 - 

inline  int  Line (//<===================================== LINE -GRID  INTERSECTION 

int  m,int  n,//< - NUMBERS  OF  GRID  INDICES  (m  AND  n  SHOULD  BE  AT  LEAST  2) 

const  double  L[6],//<--A  SCALED  LINE  SEGMENT  {L0SX, L0SY, L0Z, L1SX, L1SY, LIZ} 

int*A, int*B){//< - STORAGE  FOR  CELL  INDICES  (FOR  EACH,  SIZE  >=  m+n-3) 

int  i=SafeIndex(m,*L), j=SafeIndex(n, L [ 1] ),I=SafeIndex(m, L[3] ), 
3=SafeIndex(n,L[4]),di=i<I?l: -l,dj=j<3?l:-l,bi=di>0?l:0,bj=dj>0?l:0,k=0; 
double  M=0-j?(L[3]-*L)/(L[4]-L[l]) :0,N=I-i?l/M:0,P=*L-L[l]*M,Q=-P*N; 
if (abs(I-i)>abs(D-j))for(int  it;dj*(D-j)>=0; j+=dj,i-=di)//dj*  selects  <  or  > 
for(it=j-D?SafeIndex(m, ( j+bj)*M+P) :I;di*(i-it)<=0;i+=di)A[k]=i,B[k++]=j; 

else  for(int  jt;di*(l-i)>=0;i+=di,  j  -  =d  j ) // . di*  selects  <  or  > 

for( jt=i-I?SafeIndex(n, (i+bi)*N+Q) :D;dj*(j-jt)<=0; j+=dj)A[k]=i,B[k++]=j; 

return  k; // . the  number  of  elements  to  check  in  A  and  B 

}// - YAGENAUT@GMAIL .  COM - - LAST~UPDATED~16AUG2014 - 

templatecclass  T>void  Gradient (//< - SURFACE  GRADIENT 

const  T&Z, / / < - POINTER  TO  Z  WHERE  Z[i]  [  j]=Z(X[i] ,Y[  j] ) 

int  i,int  j,//< - GRID  INDICES  (0<=i<m-l  &  0<=j<n-l) 

double  sx, double  sy,//< - A  GRID  LOCATION  IN  SCALED  COORDINATES 

double&delsx, double&delsy){//< - SCALED  GRADIENT  COMPONENTS  (CALCULATED) 

delsx=(sy-j-l)*(Z[i][ j] -Z[i+1] [ j] )-(sy-j)*(Z[i] [ j+l]-Z[i+l] [ j+1] ); 
delsy=(sx-i-l)*(Z[i] [ j] -Z[i] [ j+1] )-(sx-i)*(Z[i+l] [ j] -Z[i+1] [ j+1] ); 

}// - YAGENAUT@GMAIL .  COM - LAST~UPDATED~16AUG2014 - 


#endif 


Topography  Example  -  The  Far  Side  of  the  Moon 

#include  "y_bmp.h"// . yBmp 

#include  "y_bilinear .  h"// . yBilinear, <cmath>{f abs ( ) ,  sqrt ( ) } 

int  main( ){//<==============CREATE  VARIOUS  IMAGES  BASED  ON  MOON  TOPOGRAPHIC  DATA 

// - READ  ORIGINAL  IMAGE  AND  CAPTURE  SCALE - 

unsigned  char *I=y Bmp: :ReadBmpFile("604359main_WAC_CSHADE_O000N1800_1000. bmp"); 

double*X=new  double[582];// . scale  values  [0,1] 

unsigned  char**C=new  unsigned  char*[582]j//. scale  colors  ({B,G,R}  and  [0,255]) 
for (int  k=0;k<582;++k)X[k]=k/581. ,C[k]=yBmp: :GetPixel(I, 1060,k+143); 

// - CAPTURE  AND  INTERPRET  A  PORTION  OF  THE  ORIGINAL  IMAGE - 

double  Z[600]  [300]; // . elevation  values 

for(int  p=0;p<600;++p)for(int  q=0;q<300;++q){ 
unsigned  char*c=yBmp: :GetPixel(I, p+200, q+400); 
double  d,D=-l; 
for(int  k=0;k<582;++k){ 

d=(C[k] [0]-c[0] )*(C[k] [0]-c[0])+(C[k] [l]-c[l])*(C[k][l]-c[l]) 
+(C[k][2]-c[2])*(C[k][2]-c[2]); 
if (d<D | |D<0)D=d,Z[p][q]=X[k];}} 

// - CREATE  A  NEW  (MONOCHROME)  IMAGE  BASED  ON  ELEVATION  VALUES - 

unsigned  char*I2=yBmp: :NewImage(600,300,255); 
for(int  p=0;p<600;++p)for(int  q=0;q<300;++q){ 
unsigned  char*c=yBmp: :GetPixel(I2,p,q); 
c [0]=c[l]=c[2]=int(Z[p] [q]*255) ; } 
yBmp : : WriteBmpFile( "moon_topography . bmp" , 12) ; 

// - CREATE  A  SURFACE -OCCULTATION  IMAGE - 


double  L[6]={0,0, 5};// . vantage  point 

int*A=new  int[600+300-3], *B=new  int [600+300-3]; 
for(int  p=0;p<600;++p)for(int  q=0;q<300;++q){ 

L[3]=p,  L[4]=q,  L[5]=Z[p]  [q]+. 00000001;// . viewable  point? 

int  K=yBilinear: : Line(600,300,L,A,B),Q=0; 
double  x,y,z; 


for (int  k=0;k<K&&Q!=4;++k)Q=yBilinear: :CellIntersect(Z,A[k],B[k],L,x,y,z); 
if (Q==4){ 

unsigned  char*c=yBmp: :GetPixel(I2,p,q); 


c[0]=c[l]=c[2]=0; }} 

yBmp : : WriteBmpFile( "moonoccultation . bmp" , 12) ; 

// - CREATE  SURFACE-GRADIENT  IMAGE - 

for(int  p=0;p<600;++p)for(int  q=0;q<300;++q){ 

int  i=yBilinear: :SafeIndex(600,p), j=yBilinear: :SafeIndex(300,q); 
double  delx,dely;/*<-*/yBilinear : :Gradient(Z,i, j,p,q,delx,dely); 
Rainbow(yBmp: :GetPixel(I2,p,q), sqrt(delx*delx+dely*dely) ,0, .4); } 
yBmp : : WriteBmpFile( "moongradient . bmp" , 12) ; 

}// - YAGENAUT@GMAIL.COM - LAST-UPDATED-16AUG2014 - 


604359main_WAC_CSHADE_O000N1800_1000 . bmp 


moon_topography . bmp 


moon_gradient . bmp 
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